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Meso-scale simulation of structure forma- 
tion for AB-dimers in solution W monomers 
was performed by dissipative particle dy- 
namics (DPD) algorithm. As a simu- 
lation model, modified Jury Model was 
adopted [Jury, S. et al. "Simulation of am- 
phiphilic mesophases using dissipative par- 
ticle dynamics," Phys. Chem. Chem. Phys. 
1 (1999) 2051 2056], which represents me- 
chanics of self-assembly for surfactant hex- 
aethylene glycol dodecyl ether (C^Eg) and 
water(H20). The same phase diagram as 
Jury's result was obtained. We also found 
that it takes a longer time to form the hexag- 
onal phase (Hi) than to form the lamellar 
phase (L a ). 

Keywords: Amphiphilic molecule, Nonionic surfactant, 
Dissipative particle dynamics simulation, Self assembly, 
Structure formation, Phase diagram 



INTRODUCTION 

Dissipative particle dynamics (DPD) simulation 
has become one of the most powerful algorithm 
for soft-matter research [1-4]. In DPD simula- 
tion, the conservative interaction potentials are 
soft-repulsive, which make simulation of long-time 
phenomena possible. By the way, DPD algo- 
rithm might be considered one of the coarse-grained 
method of molecular dynamics(MD) simulation. 
Some information of interaction potential between 
particles are neglected and simplified. We, there- 
fore, must pick up the dominant interaction poten- 
tial for the mesoscopic structure formation. Be- 
cause we don't have enough experimental data for 
the interaction potentials, it becomes a difficult 
problem of DPD simulation how we define inter- 
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action parameters. In 1999, Jury et al. suc- 
ceeded, by an empirical method, the DPD simu- 
lation of the smectic mesophase for a simple am- 
phiphilic molecule system with water solvent [5]. 
They showed that their minimal model (we call it 
Jury model) , which is composed of rigid AB dimcrs 
in solution W monomers, is very proper to present 
phase diagram of surfactant hexaethylene glycol 
dodecyl ether (Ci2E 6 ) and water (H 2 0) [5,6]. In 
this paper, we reveal processes of self-organization 
of one smectic mesophase by modified Jury model 
which has such a difference from the original Jury 
model that attractive harmonic oscillator potential 
is used as a covalent bond potential between A and 
B particles in the same molecule. Details of sim- 
ulation method are explained in the next section. 
In the last section, we will show such a interesting 
property of the system that it take a longer time 
to form the hexagonal phase (Hi) than the lamel- 
lar phase (L a ) which is more symmetric than Hi 
phase. It is considered intuitively that it takes a 
longer time to form the higher ordered phase than 
to form the lower ordered phase. However, this sim- 
ulation shows that it takes a longer time to form 
the lower ordered phase than to form the higher or- 
dered phase contrary to the intuitive consideration. 



SIMULATION METHOD 

DPD Algorithm 

First of all, we express the DPD model and algo- 
rithm [2,5]. According to ordinarily DPD model, all 
atoms are coarse-grained to particles whose mass 
are the same one. We define the total number of 
particles as N. The position and velocity vectors of 
particle i, (i = 1, • • • , N), are indicated by Vi and 
Vi, respectively. The particle i moves according to 
the following equation of motion, where all physical 
quantities are made dimensionless to handle easily 



in actual simulation. 

dt 
dVi 

dt 



Vi, 

N 

E 



ij, 



(1) 



(2) 



where a particle i interacts with the other particle 
j by Fij is the total force 



a ij 
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Table 1: The table of the coefficients aij, which depend 
on kinds of particles i and j; W is a "water" particle, A is a 
"hydrophilic" particle and B is a "hydrophobic" one. 



(3) 



In Eq. |31 F ^ is a conservative force deriving 
from a potential exerted on particle i by the 
j— particle, Fy and F^ are the dissipative and 
random forces between particle i and j, respec- 
tively Furthermore, neighboring particles on the 
same amphiphilic molecule arc bound by the bond- 
stretching force Ffj . 

The conservative force F has the following 
form: 



-V 4 <; 



(4) 



where Vj = d/dri. For computational conve- 
nience, we adopted that the cut-off length as the 
unit of length. It is assumed that the conserva- 
tive force F are truncated at this radius. Follow- 
ing this assumption, the two-point potential <pij in 
Eq. U|is defined as follows: 

1 2 

4>ij = 4>{n 3 ) = 2 a u {nj - 1) H(l - r tj ), (5) 



rj — n. We also define the 



where r i3 = |ry|;r 

unit vector n^ = rij/rij between particles i and j. 

The Heaviside step function H in Eq. is defined 

by 

{0 for x < 0, 
\ at x = 0, (6) 

1 for x > 0. 

Espahol and Warren proposed that the following 
simple form of the random and dissipative forces, 
as follows [7]: 



F 

F 



aui K {rij)ni 



Ci 



(7) 



^D^jJK-n.j)^-, (8) 



where Qj is a Gaussian random valuable with zero 
mean and unit variance, chosen independently for 



each pair (i,j) of interacting particles at each time- 
step and Qj — (ji ■ The strength of the dissipative 
and random forces is determined by the dimension- 
less parameter a and 7, respectively. The parame- 
ter At is a dimensionless time-interval of integrat- 
ing the equation of motion. 

Now we consider the fluctuation-dissipative theo- 
rem of DPD method. The time evolution of the dis- 
tribution function of the DPD system is governed 
by Fokker-Planck equation [7] . The system evolves 
to the same steady state as the Hamiltonian sys- 
tem, that is, Gibbs-Boltzmann canonical ensemble, 
if the coefficients of the diffusion and random force 
terms have the following relations: 



WD 

2 



(wr) 

2T 7 , 



(9) 
(10) 



where T is the dimensionless equilibrium- 
temperature. The forms of the weight functions 
wd and wr are not specified by the original DPD 
algorithm. We adopted the following simple form 
as the weighting functions [7] : 



WR.(r) = {cjdM} = w ( r ) 
Here the function ui is defined by [2,7]. 
ui(x) = (l-x)H(l-x). 



(11) 



(12) 



Finally, we use the following form as the bond- 
stretching force: 



B 



(13) 



where ^>b is the dimensionless bond-stretching po- 
tential energy. When particles i and j are con- 
nected in the same molecule, they interact with 
each others by the following potential energy: 



4>fj = 4> B (nj) 



1 2 



if i is connected to j. 

(14) 



Here ob is the potential energy coefficient and ro is 
the equilibrium bond length. 

Simulation Model and Parameters 

As the model, we used modified Jury model 
molecule to dimer which is composed of hydrophilic 
particle (A) and hydrophobic one (B) [5]. Water 
molecules are also modeled to coarse-grained parti- 
cles W. All mass of particles are assumed to unity. 
The number density of particle p is set to p = 6. 
The number of modeled amphiphilic molecules AB 
is shown as Nab, the number of water AT W . To- 
tal number of particle N = 2N&B + % is fixed 
to N — 10000. The simulation box is set to cubic. 
The dimensionless length of the box L is 



L = 



11.85631. 



(15) 



We use periodic boundary condition in simulation. 
The interaction coefficient ay in Eq. is given in 
Tabled 

The coefficient of bond-stretching interaction po- 
tential ae is adopted as follows: 



a B = 100. 



(16) 



We use the dimensionless time-interval of step as 
At = 0.05. 

As the initial configuration, all particle were lo- 
cated randomly. The velocity of each particle are 
distributed to satisfy Maxwell distribution with di- 
mensionless temperature T. 

The dimensionless strength of dissipative and 
random forces are 7 = 5.6250 and a = 3.3541\/T. 



total potential energy of amphiphilic molecules de- 
creases during the process from (a) to (b), more 
than from (b) to (d). (See Figd) Time depen- 
dence of the total potential energy is approximated 
to 050% (*) = 1-3389 x 10 5 - 0.081632 t, by least 
squares method for t > 200. 

For cab = 65%, the lamellar phase (L a ) was con- 
structed at i w 195 like Fig. [3] Time dependence of 
the total potential energy of amphiphilic molecules 
065% is approximated to 065% (t) = 1-2931 x 10 5 — 
0.0095747 t. (See Fig.d) Comparing 65 %(i) with 
050% (t), it is found that, the lamellar phase (L Q ) for 
cab = 65% is stable at t = 195, though structure 
for cab = 50% at t = 195 does not reach the stable 
phase, i.e., Hi phase. In other words, it takes a 
longer time to form the hexagonal phase (Hi) than 
the lamellar phase (L Q ). 

Phase Diagram 

In the previous subsection, we showed the dy- 
namics of structure formation for two cases, i.e., 
(c A b,T) = (50%, 0.5) and (65%, 0.5). Next, we 
simulated the other cases (cab,T) to obtain phase 
diagram of AB dimers in W monomers in Fig. 0J 
The phase diagram (Fig. 0J is qualitatively con- 
sistent with Jury's result [5]. The phase diagram 
(Fig. 0J also agrees with experimental result [6]. 
However, we considered temperature dependence 
of covariant bond interaction, which denotes that 
a length of AB dimers is changeable by tempera- 
ture. From this property, our modified Jury model 
is more sensitive to a change of temperature than 
original Jury model where AB dimers are rigid [5] . 



SIMULATION RESULTS 



Dynamics of Structure Formation 

We plot time evolution of total potential energy 
of amphiphilic molecules for cab = 50% and 65 %, 
at T = 0.5 in Fig. [2 

For cab = 50%, supramolecule structure grows 
to hexagonal phase around t = 900. (See FigEl ) 
Amphiphilic molecules start to make a self- 
aggregation from random configuration (FigP^a)). 
Then, local order grows like FigJ2fb). As time goes, 
the local order coalesces and the cylinder micelle 
structure grows up like FigJ^fc). At t « 900, hexag- 
onal phase (Hi) was constructed like Fig[2d). The 



1 381100 





136UO0 




! 34IMM 


1 




= 






1 ! Omim : . 


D. 




1= 




£ 


L3QIIM1 



I2S(MI» 



X ns> 



fa) ^p -1 J»9e»« - OOI&3-2/ 



c Aa = so% 




'<S) .(h) 



<■*„-«* 




-1 ^3le^0D ■ 0. OOS57 J7f 



i\ 



2M 



41M1 



M)(l 



HIM! 



di mention k j ss lime; / 



Figure 1 : Time evolution of total potential energy. Green and red dots denote total potential energy for concentration 
of AB dimes cab = 50%, and 65 %, respectively. Dimensionless temperature is kept to T = 0.5. We make snapshots at 
eight times, i.e., t =0 (a), 450 (b), 600 (c), 900 (d) for c AB = 50%, and t =0 (e), 30 (f), 75 (g), 195 (h) for c AB = 65%. 
The snapshots are drawn in Figs. 121 and I3l Using the least squares method, we fit data in the region of t > 200 to liner 
functions. 
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Figure 2: Snapshots of time evolution of self-assembly for cab = 50% at T = 0.5. As the initial configuration(a), AB 
dimers and W monomers are distributed randomly in the simulation box. Velocities of all particles arc distributed following 
Maxwell distribution. To make self-assembled structure of supramolecule clearly understandable, we pick up surface of 
the self-assembled structure, and interpolate the surface as if surface were continuous body. As time passes, hexagonal 
structure (Hi) grows like t =450 (b), 600 (c), 900 (d). 






Figure 3: Snapshots of time evolution of self-assembly for cab = 65% at T = 0.5. Visualization method, interaction 
forces, system size and other physical quantities except concentration of AB dimes are the same values as Fig. 1. As time 
passes, lamellar structure (L a ) grows like t =0 (e), 30 (f), 75 (g), 195 (h). The time to form the self- assembled structure 
is shorter than Hi phase. 
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Figure 4: Phase diagram of AB dimers in solution W 
monomers. Horizontal axis is concentration of AB dimers 
Cab- Vertical one is dimensionless temperature T. Filled 
circles denote the micelles phase (Li). Hexagonal (Hi) and 
lamellar (L a ) phases are indicated by filled diamonds and 
squares, respectively. This phase diagram is qualitatively 
consistent with the previous work by Jury et al. [5]. 



RESULTS AND DISCUSSION 

We found that it takes a longer time to form the 
hexagonal phase (Hi) than the lamellar phase (L Q ) 
where the system has higher order than Hi phase. 
This behavior is different from ordinary structure 
formation, where it takes a longer time to form the 
higher ordered phase than to form the lower or- 
dered phase. 

We would reach the conclusion that, in the pro- 
cess of structure formation in lamellar phase, the 
lamellar structure grows up directly from randomly 
distributed configuration of AB dimers. This be- 
havior is contrary to our expected scenario that 
hexagonal cylinder structure grows first, then the 
hexagonal cylinders coalesce and finally lamellar 
appears. 

Concerning the simulation model, it is easy to ex- 
tend model of amphiphilic molecules for the chain- 
like configuration. We have already studied the 
dependence of interaction parameter of phase di- 
agram [8]. 
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